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ABSTRACT 

Stimulated by recent results by iMeru &Batd (I2010allbl) . we revisit the issue of resolution 
requirements for simulating self-gravitating accretion discs with Sm oothed Particle Hydrody- 
na mics (SPH). We show that the results by iMeru & Batd (1201 Oal) are consistent with those 
of|M eru &Batd(l2010bh ifthev are both interpreted as driven by resolution effects, therefore 
implying that the resolution criterion for cooling gaseous discs is a function of the imposed 
cooling rate. We discuss two possible numerical origins of such dependence, which are both 
consistent with the limited number of available data. Our results tentatively indicate that con- 
vergence for current simulations is being reached for a number of SPH particles approaching 
10 millions (for a disc mass of order 10 per cent of the central object mass), which would 
set the critical cooling time for fragmentation at about 15il^^, roughly a factor two larger 
than previously thought. More in general, we discuss the extent to which the large number of 
recent numerical results are reliable or not. We argue that those results that pertain to the dy- 
namics associated with gravitational instabilities (such as the locality of angular momentum 
transport, and the relationship between density perturbation and induced stress) are robust, 
while those pertaining to the thermodynamics of the system (such as the determination of the 
critical cooling time for fragmentation) can be affected by poor resolution. 
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1 INTRODUCTION 

Gravitational instabilities in accretion discs have been intensively 
studied in the last decade. Several analytical models of grav- 
itationally unstable di s cs have been proposed over the years 
dBertin & Lodatolll999l: iBalbus & Papaloizoulll999l: IClarkell2009l : 
lRafikovll2009l : lRice et aLlbOldlKrumholz & BurkertlllOlOl) . How- 
ever, a large number of important properties of self-gravitating 
discs have been deter mined through the u se of n umerical simula- 
tions jGam mie 200lVLodato & Rice 2004, 2005; Rice et al. 2005 ; 
Meiia et al., ,2005.: .Bolev et al., ,200fc ,Cai et al.. ,2006.: .Boss. ,2004 



20061: ICossins et al.ll2009l , 1201(1 just to mention a few). Such nu 



merical studies have investigated seve ral different aspects of self- 
gravitating disc dynamics. Some (e.g. lLod ato & Rice 2004, 2005; 
iBolev et alJbood : ICossins et alj|2009l) have considered the angu- 
lar momentum transport in non-fragmenting discs, with the aim of 
determining whether the long term disc evolution induced by the 
instability can be treated within a local framework. This is par- 
ticularly relevant for constructing local analytical models of self- 
gravitating discs dClarke 2009: Rafikov 2009). In some other cases 
'jGam miell200l'; 'Sice et al. 2005; Boss 2004; Mayer et alj |2004 
I2OOI Boss 2006; Mayer et al. 2007; Cossins et al. 2010), the fo- 



cus was on the conditions required for a self-gravitating disc to 



fragment into bound objects. This issue is extremely important in 
the context of planet formation theories, as this mechanism forms 
the basis of the so-called gravitational instability model for planet 
formation. More in general, such mechanism might be responsible 
for the f ormation of stellar clusters in galactic c enters, including 
our own ( lNavakshinir2006l : lNavakshin et alj2067b . The mechanism 
leading to fragmentation is fairly well understood. In a cooling 
disc, thermal saturation of the instability occurs when the ampli- 
tude of the density perturbation is large enough that the associated 
shock heating provi ded by the instability balances the cooling rate 
jCossins et al.l2009h . If the required density perturbation exceeds a 
critical value, the disc fragments. Determining the critical value for 
the density perturbation is thus equivalent to determining a critical 
cooling time below which fragmentation occurs. Such crit ical cool- 
ing tim e has been variously determined over the years: iGammiel 
( I2OOII) sets it to ?>Q,~^ (where Q, is the local angular velocity of 
the disc), based o n his two dimensional shearing sheet simulations, 
iRice et all l l2005h set it to GQ,~^ for a specific heat ratio of 7 = 5/3, 
based on their three dimensional global SPH simulations, and note 
a dependency of the critical cooling time on 7. Other studies have 
found other dependencies either on the specific form of the cool- 
ing rate Johnson & Gammi&,2003. : ,Cossins et al.„2010i) or on the 
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thermal history of the disc dClarke et alj|2007l) . It is convenient to 
express the cooling time in units of the local dynamical time in the 
disc and introduce a "cooling parameter" /3 = Q,tcooi- 

Recently, some new results have cast doubts on the reliability 
of the determination of such critical cooling time, and in particular 
on those obtained by SPH simulations. In partic ular, we consider 
here two new results: first. iMeru & Bate! ( 1201 Oal) have run a num- 
ber of simulations with different slopes of the disc surface density 
profile, different disc masses, and have considered the location at 
which the fragment appeared in their fragmenting runs. Their con- 
clusion is that the critical value of /3 at a given radius i? is a func- 
tion of the "local disc mass", m{R) = T,{R)R^ /Mi,, where E is 
the surface density and A/* is the mass of the central object. In 
particular, they find that the critical cooling parameter for f r agmen- 
tation increases roughly as m ' 2. Secondlv. lMeru & Bate! f2010bh 
have noted that, for a given disc set up (in terms of total disc mass, 
surface density profile, etc.) the critical /? for fragmentation is an 
increasing function of the total number of SPH particles used, and 
— more importantly — they find no hint of convergence. 

In this paper, we reconsider the issue of resolution require- 
ments in order to simulate fragmentation in a cooling gaseous discs. 
We introduce a new resolution requirement that relates the number 
of SPH particles used to the externally i mposed cooling rate. W e 
then show that the whole set of results by iMeru & Bate ( l2010al lbl) 



can be naturally explaine d as a consequence o f failing to satisfy 
this condition. Since the iMeru & Bate! ( l20T0allbh simulations are 
among the highest resolution SPH simulations of fragmenting discs 
to date, this failure is shared also by all previous SPH simulations. 
We will then discuss which of the previous results are reliable and 
which should be considered with more care. 



2 RESOLUTION REQUIREMENTS FOR 
SELF-GRAVITATING COOLING DISCS 

Several papers in the past have investigated the resolution re- 
quirements for resolving gas fragmentation in numerical simula- 
tions. In particular. Bate & Burkert ( 1997) require that the mass 
of SPH particles within a smoothing sphere be less that the Jeans 
mass. A similar condi tion, applicable also to grid-based simulations 
( lTrueloveetaLlll997h . is that the minimum resolvable length (i.e. 
the mesh size or the smoothing length h for grid-based and SPH 
simulations, respectively) be less than the Jeans length, although 
note that failing to satisfy the Truelove condition in grid-based sim- 
ulations generally induces spurious fragmentation, while in SPH 
poor resolutio n may well lead to suppression of fragmentation. 
lNelsonl (l2006'l lists three different conditions: the first is that adap- 
tive gravitational softening is used (this is actually usually done in 
most modem SPH simul ations, including those disc ussed here), the 
second is a variant of the lBate & Burkera i ll 9971) and lTruelove et al.l 
( Il997h conditions, and the third is that the smoothing length h be 
less than the dis c thickness H. Howeve r, als o this last cond i tion is 
equivalent to the lBate & BurkertI ( 1 1 9971) and lXruelove et all ( 1 1 9971) 
condition for a marginally stable disc, for which the Jeans length is 
actually equal to the disc thickness. All these conditions can there- 
fore be cast in the simple requirement: 



h_ 
H 



< 1. 



(1) 



The natural measure of resolution for an SPH simulation of 
an accretion disc is thus the quantity h/H. In order to discuss the 
outcome of gravitationally unstable discs we are interested in eval- 



uating this quantity at the onset of the instability, that is in a con- 
dition where the disc is marginally stable. Let us then estimate the 
ratio h/H foi a marginally stable disc. We first recall than marginal 
stability for a Keplerian disc implies that 
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1, 



where Cs is the sound speed. This is equivalent to 



H 



ttQ- 



Mi, 



nQm{R), 



where m{R) = ER^/M^,. 

The SPH smoothing length is set by the condition 

7,3 3 

ph — T) nip, 



(2) 



(3) 



(4) 



where p is the local SPH density, rjisa numerical parameter gener- 
ally set to 1.2 and rUp = M^i^c/N is the SPH particle mass, where 
J\fdisc is the total disc mass, and TV is the total number of SPH par- 
ticles. With this prescription, a smoothing sphere contains roughly 
60 SPH "neighbours". We evaluate the smoothing length at the disc 
midplane. The midplane density pc is related to the surface density 
by pc — E/2_ff, which, inserted in Eq. (|4j, and using Eq. l[3) and 
the definition of m{R), gives: 
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where q — Mdisc/A^*- Actually, if the disc is very poorly resolved, 
that is if h/H > 1, the SPH method underestimates the midplane 
density by a factor H/h, and we have pc ~ p{h/H). This corre- 
sponds to an equivalent 2D version of equation Q: 



E/i^ ~ 2rf'mp. 



(6) 



In this poorly resolved case, one readily gets that h/H is given by 
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It is also useful to evaluate h/H at the outer edge of the 
disc. Rout- If the surface density profile is a power-law with re- 
spect to R with index p < 2, we have that q = Mdisc/M* = 
27rm(i?out)/ (2 — p), from which we get (for h/H < 1) 
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^2.5 X 105 

where the last equality holds for p = 1. The equivalent of (|8) when 
h/H > 1 is 

1/2 
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TV 
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(9) 

^2.5 X 105 ' ^ ' 

Note however that equation ^ is of little practical use, since in 
order to have h/H > 1 at the disc outer edge, the number of SPH 
particles needs to be very small, TV < 6000. 
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We can now reconsider the results of lMeru & Bat3 feoiodlbh 
in tiie light o f Equations ([5]> and ([8) (or ^ and In fact, on 
the one hand* Meru & BatabOlOal) find that, for a given number 
of particles, there is a relation between the imposed value of /? and 
the quantity m(i?), evaluated at the location of the first fragment 
appeari ng in their simulations. On the other hand, iMeru &Batel 
l l2010bh demonstrate that, for a given disc setup, the critical value 
of P for fragmentation is an increasing function of A'^. Both re- 
sults can be thus explained if (a) fragmentation is suppressed for 
low resolution and (b) the resolution requirement, as measured by 
the paramet er h/H, is a function of the imposed cooling time. In 
this picture, iMeru & Batd ( 1201 Oal) thus determine the location at 
which, for a g iven A'^, the resolution requirement is satisfied, while 
|m eru & Batel fcoiOh) determine the minimum A'^ at which, for a 
given setup, the requirement is satisfied. 

In fa ct, it would be puzzling to interpret the results of 
iMeru & B ate (2010a) otherwise. By construction, such simulations 
are self-similar and, as long as global effects do not play a role 
(which is the case when m{R) < 0.1), there should be no preferred 
location in the disc. The only possible mechanism that breaks the 
self-similarity is in fact resolution. 

Stated otherwise, we can say that in order to detect fragmen- 
tation in a cooling disc, the resolution requirement must depend on 
/3, in such a way that discs with longer cooling times require a bet- 
ter resolution (that is, a smaller h/H). In this case, fragmentation 
would then first occur at the minimum radius for which such resolu- 
tion requirement is satisfied, since the dynamical and cooling times 
are shorter than at larger radius. We can express such a condition in 
the form 



/3 < I3rcs= Prcsih/H) 



(10) 



What is the actual functional relation between /?rcs and h/ HI 
We explore this by plot ting in Fig. [Ti t he valu es of jS in the frag- 
menting simulations of iMeru & Bate! ( 1201 Oal) versus H/h at the 
point of fragmentation, and assuming that the correlation in these 
two quantities is purely driven by resolution, that is that /? = Pias- 
The quantity H/h is computed using equation (|5j or l[7]l, and us- 
ing the data provided by Meru & Bate (2010a) in their table 2. We 
also indicate with a leftward pointing error bar the cases where our 
estimate of m{R) (which is based on the initial disc surface den- 
sity profile) is an over-estimate: in these cases fragmentation occurs 
close to the disc inner edge and the surface density has been signif- 
icantly depleted by boundary effects by the stage of fragmentation. 

There are evidently a variety of functional forms that will pass 
through the data - we are hampered both by the fact that there is 
considerable scatter and by the fact that we have no simulation 
points with H/h > 3. One possibility is that there is a linear rela- 
tionship between (5 and H/h (see Section 3 for a discussion of this 
form in terms of the role of artificial viscosity). We have fitted the 
data with a least squares fit (shown with the solid line in Fig.[Tll and 
have obtained: 




H/h(R 



frag/ 



Figure 1. Correlation between linear resolution at the location of fragmen- 
tation, as measured by h/H, and the imposed cooling time, as measured by 
the parameter /3. Data points are taken from table 2 of Meru & Bate ( 2010ii|). 
A leftward pointing error bar indicates cases where fragmentation occurred 
near the disc inner boundary. The solid line shows a simple least squares 
fit to the data, assuming a linear relation (Eq. (TT)). The best fit model in 
this case is /3 = /3res = 2{H/h). The dashed line shows a least squares fit 
to the data assuming Eq. )121 . The best fit parameters are /3o = 14.7 and 
a = 1.77. 



expect that at high H/h the data would start deviating from the 
relation /3 = /3res and will saturate at a finite value of /3. 

Alternatively, the dashed line shows a fit to the data of the 

form 



/3res — 



{1 + ah/HY 



(12) 
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res ~ 2 I ^ 



(11) 



(see Section 3 for a discussion of the numerical effects that could 
give rise to such a relation). The important difference between this 
functional form and that given by equation dl It is that it implies 
that /3rcs converges to a finite value (/?o) at high H/h. The best fit 
values of the parameters are /3o = 14.7 and a = 1.77. We have 
excluded from the fit the data points that correspond to completely 
unresolved situations, that is h/H > 1. 

Evidently, either of the lines shown in Figure 1 provide an 
adequate fit to the data over the limited dynamic range of H/h 
available^ 

W e now turn our attention to the results of iMeru & Batd 
( l2010h) . who had found an increase in the critical j3 for fragmen- 
tion at any radius with A*'. These authors consider the case p = 1 
and q — 0.1. In the case where we are simply interested in frag- 



It is worth stressing that if the trend /3 = Acs, given by equa- 
tion dllb . extended to indefinitely high H/h then it would mean 
that there would be no physical limit to the cooling rate required 
for fragmentation. Instead it would imply that - however slow the 
cooling rate in the disc - it would always fragment if modeled at 
high enough resolution. In contrast, if there is indeed a physical 
limit to the cooling rate required for fragmentation then one would 



We also note that the data is adequately fitted by the suggested parameter- 
isation of Meru & Bate (2010a), in which /3res (in scaling with m(R)^'^) 
effectively scales with {H/h)^'^ . It is also worth noting the slightly differ- 
ent dependence of our coordinate H/ h on stellar mass compai'ed with that 
of MeiTj & Bate: we note that all disc quantities depend on the central object 
mass only via the Keplerian angular velocity, f2 and that consequently /3rcs 
mu.st depend on m{R) and q in the combination m(R)/q^^^ . 
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mentation and not in the exact location within the disc where this 
happens, we only need to make sure that the resolution criterion is 
satisfied at the disc outer edge, that is 



2/3 



iV 



2.5 X 105 



1/3 



(13) 



where the last approximate equality holds in the case that we adopt 
the linear fit (equation 1 11 1), and where we have used equations {8} 
and Onil. 

The data points of iMeru & Bat3 j2010bh are plotted in Fig.O 
along with our prediction for /3rcs based on Eq. ( I13t (solid line) and 
also the corresponding relationship based on equation M2\ (dashed 
lin e). We obtain t he thres hold for fragmentation /3frag from the data 
of lMeru&Bat3 ( l2010bl) by averaging the highest value of their 
fragmenting runs with the lowest non fragmenting run. For this pur- 
pose, we consider as non fragmenting those runs that were classi- 
fied as 'borderline' bv .Meru & Bate (2010b)^; 

We immediately see that the results of iMeru & Bate! ( 1201 Oal) 
and lMeru & Bate! j2010bh are mutually consistent if both sets of re- 
sults are explained in terms of resolution effects (that is, if /3frag = 
Pics), rcgardlcss of the exact parameterisation that is used to de- 
scribe how the cooling rate for fragmentation depends on H/h. 
We also note that, fo r the highest number of particles adopted by 
lMeru&Bat3 ( l2010bl) , the results appear to deviate from the A'^ ' 
(solid) line implied by the linear fit (equation II It. In these cases, 
we thus have non-fragmenting simulations that are well resolved 
according to this criterion. 

The simulation results in Figure 2 are also broadly consistent 
with the dashed line (constructed using equation l ll2t ) which is 
a parameterisation that implies a converged value of the required 
cooling rate for simulations that are sufficiently well resolved. 

Figure 2, and especially the deviation from the A''^'^'^ depen- 
dence at high A'^, therefore contains tentative evidence that conver- 
gence may be being reached at the largest number of particles, and 
that the 'true' value of the threshold cooling time for fragmentation 
would therefore lie between /3 « 10 — 15, that is roughly a factor 
of two larger than previously thought. 




Figure 2. Correlation between the critical cooling time for fragmentation 
and the number of SPH particles used, for the simulations described in 
Meru & Bate 1 2010b). The threshold is obtained averaging the highest /3 
for fragmenting simulations and the lowest /3 for non-fragmenting ones. 
The solid and dashed lines show equations (TT) and jl2t . respectively, eval- 
uated at the disc outer edge, and where the free parameters are set by the fit 
shown in Fig.[T] 



sure that such external heating sources, and in particular numerical 

heating sources, are negligible. 

It is well known (Pringle 1 98lh that the condition of thermal 
equilibrium for a viscous disc can be easily e xpressed in terms of 
a rela tion between the viscosity parameter a jShakura & SunvaevI 
Il973h and the cooling parameter /?: 



3 PHYSICAL ORIGIN OF THE RESOLUTION 
CONDITION 

In the previous section, we have explained the available simulation 
data in terms of two alternative parameterisations of how the cool- 
ing rate required for fragmentation may depend on the ratio of the 
smoothing length to the disc scale height. We now turn our atten- 
tion to possible physical explanations of such requirements. 

The first possibility is related to artificial viscosity. Indeed, 
the artificial visco sity term in SPH is linearly proportional to h/H 
( lMonagharJl992l) , and its effect might in principle alter the thermal 
balance of the disc. 

As discussed above, for cooling discs the saturation of the 
gravitational instability is due to the balance of internal heating 
due to shocks induced by the instability itself and the externally 
imposed cooling. Clearly, if any other form of external heating is 
present, the saturation amplitude of the instability will be decreased 
and fragmentation might therefore be inhibited even for cooling pa- 
rameters that would lead to fragmentation in the absence of exter- 
nal heating. For actual astrophysically relevant discs, such external 
heating sources might well be an important physical ingredient to 
be taken into account. On the other hand, in numerical simulations, 
where one wishes to determine the critical cooling rate, one has be 



97(7 



(14) 



For SPH simulations, the most significant source of numeri- 
cal heatin g is provided by artificial viscosity. In particular, it can be 
show n jArtymowicz & Lubowll996lJMurravll996l : lLodato & Pried 
I2OIOI) that the linear term in the standard artificial viscosity imple- 
mentation of SPH leads to an equivalent Shakura-Sunyaev param- 
eter Oart given by 



1 h 

Qart — T7:Q!SPH-77, 

iU H 



(15) 



where qsph is an SPH input parameter (usually set to 1, but re- 
duced to 0.1 in the simulations discussed here). Actually, the above 
relation holds only in the continuum limit and if artificial viscosity 
is applied both in regions of converging and diverging flow. Usu- 
ally, a switch enabling artificial viscosity o nly for converg ent flows 
is used, and sometimes additional switches ( lBalsarall99a) are used, 
so that the above relation is actually an overestimate of the resulting 
viscosity. 

When ftart becomes comparable to the value implied by Eq. 
( 114b . artificial heating offsets completely the imposed cooling, and 
no gravitational perturbations will be able to grow. This occurs 
when 
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where the last equality refers to the case of 7 = 5/3 and 
agPH = 0.1. The resolution criterion proposed in the previous 
section (equations [TO] and IIU thus corresponds to the case where 
artificial viscosity provides 5 per cent of the heating required for 
thermal equilibrium. 

We draw attention to two aspects of this result. Firstly, it would 
imply - rather surprisingly - that artificial viscosity effects can pre- 
vent fragmentation even when they contribute a tiny fraction (5%) 
of the total thermal energy balance. Secondly, it should be stressed 
that this condition is a necessary (but not necessarily a sufficient) 
criterion for fragmentation. As noted in Section 2, there is tenta- 
tive evidence from Figure 2 that the points at the highest fall 
below the solid line - in other words we have simulations that are 
not fragmenting even though they are in the regime where the arti- 
ficial viscosity contribution is less than 5%. We interpret this result 
as evidence that there is indeed a physical mechanism preventing 
fragmentation for small cooling rates. 

In addition to the thermal effects described above, artificial 
viscosity might also dynamically stabilize the disc. The specific 
effect of viscosity in this regard depends on how viscosity scales 
with surface density and mi ght also lead to a secular instability 
jSchmit & TschamuteJfTggl) . While this is beyond the scope of 
the present paper, such effects should be further investigated. 

Since we find it rather surprising that artifical viscosity should 
suppress fragmentation when it contributes such a minor compo- 
nent to the disc's thermal balance, we also explore the hypothesis 
that the effect of under resolution is one of artificially smoothing 
density enhancements over the SPH smoothing length h. In poorly 
resolved simulations, where h is not much smaller than the length 
scale of density peaks (A), then this will suppress the peak ampli- 
tude of the resulting density fluctuations by a factor of 1 + /i/A (as- 
suming that gravitational in stabilities generate predominantly lin- 
ear structures). According to lCossins et the gravitational 
heating rate associated with modes of r.m.s. fractional amplitude 
AE/E is proportional to (AE/E)^ and therefore in thermal equi- 
librium this quantity scales with the cooling rate (oc /3~^). If we 
assume that in a perfectly resolved calculation the peak amplitude 
would scale with the r.m.s. amplitude (but be degraded by a factor 
1 + h/X in the case of finite h) and if we furthermore associate 
fragmentation with the peak fluctuations achieving a critical value 
of AE/E (of order unity) then it follows that the cooling rate re- 
quired for fragmentation is increased in the case of simulations of 
finite h. Specifically we have 

/3rcs = + (17) 

where Po is the value required for fragmentation in the case of a 
well resolved simulation (i.e. small h). In estimating A we note that 
the Jeans length in a disc with Toomre Q parameter close to unity 
is just the vertical scale height H; it therefore seems logical that A 
should scale with H and this motivates the form of equation | |I2I >. 

We note that both our resolution criteria do not abandon the 
notion that (in a well resolved simulation) there is a critical cooling 
rate associated with fragmentation but take into account that the 
fact that more vigorous cooling (lower /3) is required in the case of 
poorly resolved simulations. 



4 DISCUSSION AND CONCLUSIONS 

In this paper we have shown that the results of iMeru & Batd 
( l2010al) (concerning the location at which fragmentation occurs 
in SPH simulations of cooling self-gravitating discs) are quanti- 
tatively consistent (if in terpreted as being dri ven by resolution ef- 
fects) with the results of lMeru & Batd ( 1201 Obf) in which it is shown 
that the critical cooling rate for fragmentation depends on the num- 
ber of particles . Indeed the results of each paper imply the results 
of the other and both imply that the measured cooling rate for frag- 
mentation is a function of h/H (where h is the SPH smoothing 
length and H is the disc scale height). 

It is difficult to assign a precise functional form for the way 
that the critical cooling rate depends on resolution, given the scat- 
ter in the simulation data and the relatively small range m H/h 
for which simulations are currently available. We have explored 
two possibilities for effects that may explain the simulation results. 
In one case we explore the possibility that a necessary criterion 
for fragmentation to be detected is that the artificial viscosity con- 
tributes less than a certain fraction of the thermal energy input (and 
associated angular momentum transport) in the disc. We find that 
the results at low resolution can be explained in these terms but 
that the requirement on the contribution from artifical viscosity 
is suprisingly stringent (i.e. fragmentation occurs only where this 
contributes less than 5% of the energy input to the disc). Alterna- 
tively, we consider the possibility that finite resolution just smooths 
out density peaks so that fluctuations that would collapse if well re- 
solved do not achieve a critical amplitude at finite h. 

Irrespective of how one explains how resolution affects the 
cooling req uirements for fragm entation, an important issue - as 
stressed by iMeru & Batd |2010b) - is whether there is indeed a 
convergence in the required cooling rate at high resolution. In other 
words, is there a level of resolution above which one recovers the 
result that has been assumed hitherto - i.e. that there is a physical 
condition on the cooling rate required for fragmentation? We evi- 
dently cannot answer this question with calculations that have not 
attained convergence (see Figure 2 where /3frag rises mildly with 
TV even at the highest A'^ values studied). However, Figure 2 con- 
tains hints of approaching convergence. One way to see this is to 
compare the simulation data with the solid line which corresponds 
to requirement of a 5% contribution from aritifical viscosity. The 
fact that the right hand points lie below this solid line (where /3frag 
scales with N^^"'') is evidence that this condition may be necessary 
but not sufficient for fragmentation. Similarly, the simulation data 
is consistent with the dashed line in which the actual /? required for 
fragmentation converges to a value of /3o ~ 14.7 at high resolution. 

Clearly the question of ultimate convergence will only be set- 
tled by further simulations at higher A'^. If convergence is not ulti- 
mately attained then it would involve a radical re-think of all our 
understanding of gravitationally unstable discs since it would im- 
ply that any self-gravitating disc (however slowly cooling) should 
in reality fragment. This would be surprising because, in the case 
of very slow cooling, the r.m.s. mode amplitude required to achieve 
thermal balance should be very low and it would be unexpected 
(though possible in principle) that such a disc nevertheless exhib- 
ited locally non-linear density fluctuations giving rise to collapsing 
fragments. 

The more conservative conclusion (if convergence is ulti- 
mately attained) is that we have simply under-estimated the critical 
/3 value hitherto. This would involve some quantitative adjustment 
of previous conclusions. For example, if the critical cooling time 
for fragmentation is roughly a factor 2 higher than had been previ- 
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ously thought then this affects the locatio n at which giant pla nets 
can form through gravitational instability jMeru & Batel2010bll : in 
the regime of ice cooling appropriate to the outer parts of protostel- 
lar discs, this change would bri ng in the min i mum radius for frag - 
mentation by a factor ~ 2^/^ ( IClarkell2009l : ICossins et alJIlOlOh . 
Although this change is quite modest, it may be important to the 
discuss ion of whether re cently imaged planets (e.g., around H R 
8799, Marois et al.' 2008, around /3 Pic, iLagrange et all |2009| or 
around Fomahault, Kalas et al. 2008) could have been formed by 
gravitational instability. Note, however, that — for a given critical 
cooling time — the uncertainty in the determination of the frag- 
mentation radius due to uncertainties in the relevant opacities, in 
the detailed vertical structure of the disc and due to the effects of 
magnetic fields are probably larger than the modest change dis- 
cussed above. 

More in general, we might ask which ones of the properties 
of non fragmenting self-gravitating accretion discs are going to be 
significantly affected by poor resolution, in the sense expressed by 
Eq. l llOt . Probably the most important property, apart from frag- 
mentation, is related to the stress induced by the instability, and in 
particular its magnitude and the locality of the associated dissipa- 
tion (which are in turn related to the power spectrum of the modes 
excited by the instability). 

As mentioned in the Introduction , sev eral studies 

jLodato & Ri^ |2004 l2005l : ICossins et all l2009h have con- 
sidered the issue of locality of the transport induced by the 
instability. Here the question is whether the typical wavelength of 
the disturbances is a signi ficant fraction of the radius R. The above 
studies (and in particular ICossins et al.ll200^ have shown that the 
spectrum of the disturbances peaks at a wavelength X ~ H and 
becomes negligible for A > 3H. In this respect, thus, the relevant 
resolution requirement is the g enerally less re strictive condition 
h/H < 1 already discussed in iNelson j2006l) . Furthermore, not 
resolving the smallest scales would only prevent the disc from 
developing small scale structures, which would not contribute to 
global transport anyway. 

Let us now discuss whether the evaluation of the stress from 
such simulations is affected by Eq. dlOt . Clearly, /or a given dis- 
turbance, regardless of the way it has been excited, a correct eval- 
uation of the resulting stress only requires that the typical wave- 
lengths of the disturbance are resolved. However, here the question 
is whether the disc response to an external boundary condition (in 
this case, to an externally imposed cooling rate) depends on res- 
olution. The condition of thermal equilibrium relates in a unique 
way the external cooling to the magnitude of the stress induced by 
the instability (eq. JI4b). under the hypothesis that the only heating 
source in the disc is provided by the instability itself and is dissi- 
pated locally. Here, we have shown that the new empirically deter- 
mined resolution condition corresponds to the case where artificial 
viscosity only provides 5 per cent of the heating necessary for ther- 
mal balance. Therefore, unless the effects of artificial viscosity are 
much larger than expected, one would not expect a change in the 
stress above the 5 per cent level. In this respect, the effect of poor 
resolution is just to suppress the development of small scale distur- 
bances, and thus reduce the peak value of the density perturbation, 
while keeping the overall average value of the stress essentially un- 
changed. 

Finally, it is worth noting that the discussi on above is fo- 
cussed on SPH simulations. Many p revious results ( lGammiel200ll : 
iMeiia et al.ll200^ ; iBolev et al.ll20ci6h . i ncluding several determina- 
tions of the fragmentatio n boundary dCammi 3 I2OOIL lMeiiaetal.1 
I2OO5I : IBossI |2004 |2006|) . have been obtained using grid-based 



codes. Artificial viscosity is also present to various degrees is grid- 
based codes, and one might thus conclude that an analogous reso- 
lution requirement should be formulated also in thus case. Whether 
also grid-based simulations do not converge in the determination of 
the critical time scale for fragmentation and the exact form of the 
resolution requirement in this case are questions that still need to 
be worked out. 
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